function [xc, lags] = xcov3(A,B)

%Thomas Stevens - Lim Lab
%opened 01/15/12
%xcov3(A,B), or xcov3(A)
%output: [normalized cross covariance, 1x3 cell of lags]
%Cross-covariance in 3D for determination of correlation length
%Default normalizes result at 0 lag to 1.
%Omission of second 3d matrix yields autocorrelation


if nargin<2 || isempty(B)
    A = A - mean(A(:));
    B=A;
else
    A = A - mean(A(:));
    B = B - mean(B(:));
end

%Use convnfft from matlab central for speed
%From doc: In 3D, this function is faster than CONVN for nA, nB > 5.
%tic
%xc = convn(A,flipall(B));
%toc
tic
disp('Convolution by FFT/iFFT...')
xc = convnfft(A,flipall(B));
toc
zero_lag = xc(floor((size(A,1)+size(B,1))/2), ...
              floor((size(A,2)+size(B,2))/2), ...
              floor((size(A,3)+size(B,3))/2));
xc = xc/zero_lag;
lags = cell(1,3);
for i=1:length(lags)
    lags{i} = floor((-size(A,i)-size(B,i))/2)+1 : 1 : ...
               floor((size(A,i)+size(B,i))/2)-1;
end

end